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Abstract In this paper we investigate two numerical schemes for the simulation of 
stochastic Volterra equations driven by space-time Levy noise of pure-jump type. 
The first one is based on truncating the small jumps of the noise, while the sec¬ 
ond one relies on series representation techniques for infinitely divisible random 
variables. Under reasonable assumptions, we prove for both methods IP- and al¬ 
most sure convergence of the approximations to the true solution of the Volterra 
equation. We give explicit convergence rates in terms of the Volterra kernel and the 
characteristics of the noise. A simulation study visualizes the most important path 
properties of the investigated processes. 
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1 Introduction 

The aim of this paper is to investigate different simulation techniques for stochastic 
Volterra equations (SVEs) of the form 



where G is a deterministic kernel function, <7 a Lipschitz coefficient and A a Levy 
basis on R + x M.' 1 of pure-jump type with no Gaussian part. In the purely temporal 
case where no space is involved and the kernel G is sufficiently regular on the diag¬ 
onal {(f;s) £ R+ x R + : t =s}, the existence and uniqueness of the solution Y to (1) 
are established for general semimartingale integrators in [14]. The space-time case 
(1) is treated in [5] for quite general Levy bases. In particular, G is allowed to be 
singular on the diagonal, which typically happens in the context of stochastic partial 
differential equations (SPDEs) where G is the Green’s function of the underlying 
differential operator. More details on the connection between SPDEs and the SVE 
(1) are presented in Sect. 2, or can be found in [2, 5, 20], 

Since in most cases there exists no explicit solution formula for the SVE (1), it 
is a natural task to develop appropriate simulation algorithms. For SPDEs driven 
by Gaussian noise, research on this topic is rather far advanced, see e.g. [7, 9, 21]. 
However, for SPDEs driven by jump noises such as non-Gaussian Levy bases, the 
related literature is considerably smaller, see [3] and the work of Hausenblas and 
coauthors [8, 10, 11]. The case <7 = 1 has been treated in [4]. The contribution of 
our paper can be summarized as follows: 

• We propose and analyze two approximation schemes for (1), each of which re¬ 
places the original noise by a truncated noise that only has finitely many atoms 
on compact subsets of K + x For the first scheme, we simply cut off all jumps 
whose size is smaller than a constant. For the second scheme, we use series rep¬ 
resentation techniques for the noise as in [17] such that the jumps to be dropped 
off are chosen randomly. Both methods have already been applied successfully 
to the simulation of Levy processes, cf. [1, 18]. 

• In the case where G originates from an SPDE, the crucial difference of our nu¬ 
merical schemes to the Euler or finite element methods in the references men¬ 
tioned before is that we do not simulate small space-time increments of the noise 
but successively the true jumps of the Levy basis, which is an easier task given 
that one usually only knows the underlying Levy measure. It is important to rec¬ 
ognize that this is only possible because the noise A is of pure-jump type, and 
contains neither a Gaussian part nor a drift. We shall point out in Sect. 6 how to 
relax this assumption. 

The remaining article is organized as follows: Section 2 gives the necessary back¬ 
ground for the SVE (1). In particular, we present sufficient conditions for the ex¬ 
istence and uniqueness of solutions, and address the connection between (1) and 
SPDEs. In Sect. 3 we construct approximations to the solution Y of (1) by truncat¬ 
ing the small jumps of the Levy basis. We prove in Thm. 1 their //-convergence. 
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and in some cases also their almost sure (a.s.) convergence to the target process Y. 
In Sect. 4 we approximate the driving Levy basis using series representation meth¬ 
ods. This leads to an algorithm that produces approximations again converging in 
the //-sense, sometimes also almost surely, to Y, see Thm. 2. In both theorems, we 
find explicit //-convergence rates that only depend on the kernel G and the charac¬ 
teristics of A. Section 5 presents a simulation study for the stochastic heat equation 
which highlights the typical path behaviour of stochastic Volterra equations. The 
final Sect. 6 compares the two simulation algorithms developed in this paper and 
discusses some further directions of the topic. 


2 Preliminaries 

We start with a summary of notations that will be employed in this paper. 

R+ the set [0,°°) of positive real numbers 

N the natural numbers {1,2,...} 

B a stochastic basis (£2,JP,F = (J^) /€ r + ,P) satisfying the usual hy¬ 

potheses of completeness and right-continuity 
£2,£2 £2 := £2 x R+ and £2 := £2 x R+ x W 1 where d £ N 

3§{FL d ) the Borel cr-field on R d 

./i} the collection of all bounded Borel sets of R + x W 1 

3* the predictable cr-field on B or the collection of all predictable processes 

£2 ->R 

3P the product 3? ® ^(W 1 ) or the collection all SP (g> <S^(R d )-measurable 

processes £2 —> R 

the collection of sets in 3? which are a subset of £2 x [0,k] x \—k,k] d 
for some k £ N 
p* p V 1 

L p the space Z/(f2, JF,P), p £ (0,°°], endowed with the topology induced 

by \\XWu, :=E[\X\p} 1 /p* 

L° the space L°(£2,&,P) of all random variables on B endowed with the 

topology of convergence in probability 

^loc the set of all Y £ 3P for which ||T(f,x)|| u> is uniformly bounded on 

[0, T] x W l for all T £ R+ (p £ (0,°°]) 

A c the complement of A within the superset it belongs to (which will be 

clear from the context) 

A—B {x— y: x £ A,y £ B} 

-A {—x: x£A} 

Leb the Lebesgue measure on R f/ (d should be clear from the context) 

|| • || the Euclidean norm on R f/ . 

C,C(T) two generic constants in R + , one dependent and one independent of T, 
whose values we do not care of and may therefore change from one 
place to the other 
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We suppose that the stochastic basis B supports a Levy basis, that is, a mapping 
A : Mb —> L° with the following properties: 

• A(0)=Oa.s. 

• For all pairwise disjoint sets C Mb with G Mb we have 

A (Qa,) = £a(A 0 inL°. (2) 

• (A(f2 x is a sequence of independent random variables if pj are 

pairwise disjoint sets in Mb- 

• For every B £ Mb, A (£2 xB) has an infinitely divisible distribution. 

• A (A) is ^-measurable when A € Mb andA C Q x [0,f] x Effort £ R + . 

• For every t £ R+, A £ Mb and f2o £ M t we have a.s. 

A (A n (£2o x (t,°o) x R rf )) = l i2o A(An(f2 x x R d )). 

Just as Levy processes are semimartingales and thus allow for an Ito integra¬ 
tion theory. Levy bases belong to the class of L°-valued cr-finite random measures. 
Therefore, it is possible to define the stochastic integral 

[ H(s,y)A(ds,dy) 

J R + xl rf 

for H £ M that are integrable with respect to A, see [6] for the details. 

Similarly to Levy processes, there exist two notions of characteristics for Levy 
bases: one going back to [15, Prop. 2.1] that is based on the Levy-Khintchine for¬ 
mula and is independent of F, and a filtration-based one that is useful for stochastic 
analysis [6, Thm. 3.2]. For the whole paper, we will assume that both notions coin¬ 
cide such that A has a canonical decomposition under the filtration F of the form 

A(dr,ck) = fi(df,dr)+A c (dr,ck) + J zl{| z |<i} (jt — v)(dr,ck,dz) 

+ J^t { \ zl>l] n(dt,dx,dz) , 

where B is a cr-finite signed Borel measure on R + x R rf , A c a Levy basis such that 
A (£2 x B) is normally distributed with mean 0 and variance C(B ) for all B £ Mb, 
and ,ti a Poisson measure on R + x M. d relative to F with intensity measure v (cf. [12, 
Def. II. 1.20]). There exists also a cr-finite Borel measure A on R. x R' / such that 

B(dt,dx) = b(t,x) A(dr,dx) , C(dr,dv) = c(t,x) A(dr,ck) and 
v(df,dr,dz) = n(t,x,dz) A(dr,dr) (3) 

with two functions b: R+ xl^-tR and c: R+ x —> R + as well as a transition 
kernel n from (R. + x R rf ,^(R + x M. d )) to such that n(t,x, •) is a Levy 

measure for each (t,x) £ R+ x 
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We have already mentioned in the introduction that we will assume 


C = 0 


(4) 


throughout the paper. For simplicity we will also make two further assumptions: 
first, that there exist fcgR and a Levy measure n such that for all ( t,x ) £ K+ x 
we have 

b(t,x)=b, 7t(t,x,-) = 7t and A(dr,dx) = d(t,x); (5) 


second, that 


A 


( 6 ) 


where +V is the collection of all symmetric Levy bases and Vq is the class of Levy 
bases with locally finite variation and no drift , defined by the property that 


J R |z|%|<t}^(dz)< 


°° , and bo := b — 


J R zl { \ z \< n n(dz) = 0. 


Furthermore, if n has a finite first moment, that is, 



l>i } 7z;(dz) < 


oo 

5 


(7) 


we define 


B 1 (dt,dx):=b l d(t,x), b\ :=b+ f zl {W>1} n{dz ), 

M(dr,dx) := A(dr,dx) —.Bi(dr,dx) = / z(/i — v)(dr.dx,dz) . 

Jr 

Next, let us summarize the most important facts regarding the SVE (1). All de¬ 
tails that are not explained can be found in [5]. First, many SPDEs of evolution type 
driven by Levy noise can be written in terms of (1), where G is the Green’s function 
of the corresponding differential operator. Most prominently, taking G being the 
heat kernel in R c/ , (1) is the so-called mild formulation of the stochastic heat equa¬ 
tion (with constant coefficients and multiplicative noise). Typically for parabolic 
equations, the heat kernel is very smooth in general but explodes on the diagonal 
t = s and x = y. In fact, it is only /;-fold integrable on |0.7'] x W 1 for p <1+2/d. In 
particular, as soon as d > 2, it is not square-integrable, and as a consequence, no so¬ 
lution to the stochastic heat equation in the form (1) will exist for Levy noises with 
non-zero Gaussian component. This is another reason for including assumption (4) 
in this paper. 

Second, let us address the existence and uniqueness problem for (1). By a so¬ 
lution to this equation we mean a predictable process Y € & such that for all 
(t,x) £ M+ x the stochastic integral on the right-hand side of (1) is well de¬ 
fined and the equation itself for each (t,x) £ [0,T] x W 1 holds a.s. We identify two 
solutions as soon as they are modifications of each other. Given a number p £ (0,2], 
the following conditions guarantee a unique solution to (1) in B/ oc by [5, Thm. 3.1]: 
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A1. Yq £ flf is independent of A. 

A2. a : R —> R is Lipschitz continuous, that is, there exists C £ R + such that 


\o(x)-G(y)\<C\x-y\, x,y £ R . (8) 

A3. G: (R + x R d ) 2 —>• R is a measurable function with G(f, -;s, •) = 0 for s > t. 

A4. A satisfies (3)-(6) and 

[ \z\ p n(dz) < °° . (9) 

Jr 

A5. If we define for (; t,x ), (s,y) £ R+ x R f/ 

G(t,x;s,y) := \G{t,x\s,y)\\ {p>hA ^ + \G(t,x\s,y)\ p , (10) 

then we have for all T £ R + 

sup / f G{t,x\s,y)d{s,y) <<*>. ( 11 ) 

(i,*)e[0,7’]xM‘ # -'° 


A6. For all e > 0 and T £ R + there exist k £ N and a partition 0 = to < • • • < 4 = T 
such that 

sup sup / / G(t,x\s,y)d(s,y) < e . (12) 

Apart from A1-A6, we will add another assumptions in this paper: 

A7. There exists a sequence (U n )neN of compact sets increasing to R f/ such that for 
all T £ R+ and compact sets K C R f/ we have, as N —> °°, 



Conditions A6 and A7 are automatically satisfied if \G(t,x\s,y)\ < g(t s,x — y) 
for some measurable function g and A5 holds with G replaced by g. For A6 see [5, 
Rem. 3.3(3)]; for A7 choose U N := {x £ R^: ||jc|| < N} such that for p > 1 


< 


sup 

(t,x)e[0,T]xK 

sup 

(f,jc)e[o,r]x/t 


/ I „ |G(f,x;s,y)| p d(s,y) 
Jo J(U N ) C 

f f ^ 

Jo J(u N y 


sup/ [ g p (s,y)d{s,y)< [ [ g p {s,y) d(j,y) 
xekJ o Jx-(u N ) c Jo Jk-(u n ) c 


—> 0 as N °° 
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by the fact that K — ( U N ) C 4- 0. A similar calculation applies to the case p £ (0,1] 
and the first term in r^(T. K). 

Example 1. We conclude this section with the stochastic heat equation in M. d , whose 
mild formulation is given by the SVE (1) with 


G(t,x-,s,y)=g(t-s,x-y), g(t,x) 


ex pHMI 2 /(4Q) 

(4nt) d l 2 




(14) 


for (t,x), (s,y) £ K.+ x W 1 . We assume that Yq and a satisfy conditions A1 and A2, 
respectively. Furthermore, we suppose that (3)—(6) are valid, and that (9) holds with 
some p £ (0,1 + 2/d). It is straightforward to show that then A3-A6 are satisfied 
with the same p. Let us estimate the rate r^(T,K) for T £ K.+, K := {||jc|| < R} with 
R £ N, and U N := {||x|| < N}. We first consider the case p < 1 or A £ U7. Since 
K — (U N ) C = (U N - R f for N> R, the calculations after A7 yield (L(•, •) denotes the 
upper incomplete gamma function and p(d) := 1 + (1 — p)d/ 2) 




eX P(-'"' 2 /(4'))/- ldrd , 


c(n U r . r (|,£(^)_(£Mi)' W 


^(d , ,, p(N — R) 2 

xT- p(d), - 

V 2 h 4T 


< C(r)exp ( — 


p{N-R) : 


4 T 


(N-R) 


d-2 


(15) 


which tends to 0 exponentially fast as N —> If p > 1 and A <f it follows from 
formula (13) that we need an extra summand for (T,K), namely (15) with p = 1. 


3 Truncation of Small Jumps 


In this section we approximate equation (1) by cutting off the small jumps of A. To 
this end, we first define for each N £N 


G N (t,x;s,y) := G(t,x\s,y)\ uN (y ), (f,x),(s,y) S M+ xl 


(16) 


where the meaning of the sets U N is explained in A7. Furthermore, we introduce 

Up* 


4-=(l \z\ p n(dz) 
\J[-e N ,e N ] 




L 


zt 


{p> 1,A^} 


?r(ck) 


(17) 
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where (e /v )/ V gpj C (0,1) satisfies e N —> 0 as IV —> Condition A4 implies that 
as N —>• oo. Next, defining truncations of the Levy basis A by 

A Ar (dr,ck):= [ zn{dt,dx,dz) , (18) 

J[-e N ,e N ] c 

our approximation scheme for the solution Y to (1) is given as: 

Y N (t,x) := Y 0 (t,x) + f [ G N (t,x;s,y)a(Y N (s,y))A N (ds,dy) (19) 
Jo JsJ 

for (t,x) £ K.+ x W 1 . Indeed, Y N can be simulated exactly because for all T £ R+ the 
truncation A N only has a finite (random) number R N (T) of jumps on [0.7) x U N , 
say at the space-time locations (t f with sizes jf. This implies that we have 

the following alternative representation of Y N (t,x) for (t,x) £ [0,7"] x W 1 : 

r n (t) 

Y N (t,x)=Y 0 (t,x) + £ G(t,x-T"^)o(Y N (T?^))J?t {x „ <t} . (20) 

1=1 

What remains to do is to simulate Y N (i = 1 ,... ,R N (T), iteratively, from 
which the values Y(t,x) for all other (t,x) £ [0, T] x M. d can be computed. 

The following algorithm summarizes up the simulation procedure: 

Algorithm 1. Consider a finite grid Sf that is a subset of [0,7"] x . For each step 
N proceed as follows: 

1. Draw a Poisson random variable R N (T) with intensity 

R N (T) ■= [ j uN J R t {ze[ _ eN ^]c } v(dcdMz) = TLeb (t/^)^([— , e^] c ) . 

2. For i=l,...,R N (T): 

a. Draw a pair with uniform distribution from [0,7"] x U N . 

b. Draw Jf from [—e^,e*^] 0 with distribution n/n[[—E N ,E N ] c ). 

3. For each i = 1,... ,R N (T) and (t.x) £<3 simulate Yq(t^, ^- v ) and Yo(t,x). 

4. For each i = 1,... ,R N (T) set 

j= 1 

5. For each (t.x) £l$ define Y N (t,x) via (20). 

The next theorem determines the convergence behaviour of the scheme (19) to 
the true solution Y to (1). 

Theorem 1. Grant assumptions A1-A7 under which the SVE (1) has a unique so¬ 
lution in B d oc . Then Y N as defined in (19) belongs to B d oc for all N £ N, and for all 
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T £ R + and compact sets K C there exists a constant C{T ) £ R+ independent 
ofN and K such that 

sup ||r(f,x)-F A, (f,x)|| LP <C(r)(rf(r,^) + ^ + tf). (21) 

Furthermore, ifjff =1 (T,K) + r% + r^) p < °° is fulfilled, then we also have for 

all (t,x) £ [0, T] x K that Y N (t,x) —> Y ( t,x ) a.s. as N —> °°. 

Proof It is obvious that | G n \< |G| pointwise and that we have v N <V for the third 
characteristic V :V of A N . Thus, A1-A6 are still satisfied when G and v are replaced 
by G n and v N (if A £ y, also A N £ y). So Y N as a solution to (1) with G N and A N 
instead of G and A belongs to B^ oc as well. Moreover, for all T £ R + there exists 
C(T) £ R+ independent of N £ N such that 

sup \\Y N (t,x)\\ LP <C(T) , NeN. (22) 

(tx)6[0,7’]xR rf 

We only sketch the proof for this statement. In fact, using [5, Lem. 6.1(1)] it can be 
shown that the left-hand side of (22) satisfies an inequality of the same type as in 
Lem. 6.4(3) of the same paper. In particular, it is bounded by a constant C N (T) that 
depends on N only through | | and v N , and that this constant is only increased if 

we replace |G N | and v N by the larger |G| and v. In this way, we obtain an upper 
bound C(T) that does not depend on N. 

Next, we prove the convergence of Y N to Y as stated in (21). We have 

Y(t,x)-Y N (t,x)= [ [ [G(t,x-s,y)-G N (t 1 x;s,y)}a{Y(s,y))A{ds 1 dy) 

Jo Js. d 

+ [ [ G N (t,x;s,y)[o{Y(s,y))-o(Y N (s,y))]A(ds,dy) 

Jo Jw> 

+ f f G N ( t ’ x ’ s ,y)o(Y N (s,y))(A -A^Xdv.dy) 

Jo Jw d 

=:/f(Lx)+/^(f,x)+/f(/,x) , (t,x) £ R+ x M. d . (23) 


If p > 1, we have by (8), Holder’s inequality and the Burkholder-Davis-Gundy- 
inequality 


\\l 2 {Fx)\\lp < [ [ G N (t,x-s,y)[o(Y(s,y))-a{Y N (s,y))]Bi(ds,dy) 

Jo Jw 1 

l [ d GN ^' x ’ s A)[(y(Y(^y))-(y{Y N (s,y))}M{ds,dy) 

JO J R rf 

Jo J Rd \G(t,x',s,y)\ |Bi|(d.v,dy) 


LP 


<c\ 


1 /p 


X Jo Is,d |fii|(ds,dy) 
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+ c (^J o J R J G {t, x ',s,y)\ p \\ Y (s,y)- YN (s,y)\\Lpd{s,y)^J . (24) 
If p € (0,1], we have A £ Yq by (6) and (9), and thus Jensen’s inequality gives 
\\ I ^{t,x)\\ LP =E (f f G N (t,x-,s,y)[<j{Y(s,y))-o(Y N {s,y))]zn(ds,dy,dz) 

_\J 0 JR rf 

- E / /. |G A '(fA;i,y)[c7(F(s,y))-(7(7 JV (s,y))]z| p v(ds,dy,dz) 

Jo 

< c [ [jG(t : x-,s,y)\ p \\Y(s,y)-Y N (s : y)\\ L pd{s,y) . (25) 

J 0 jR rf 

Inserting (24) and (25) back into (23), we have for v N (t,x) := \\Y(t,x) — Y N (t,x)\\u> 
v N {t : x) < C(T) ^ ^ |G(f,x;s,y)|l {p>liA ^ } (t' A, (s,y)) /, d( S ,y)^ 

+ ^J^jG(t,x;s,y)\P(^(sj)Y t d(s,y)y P j 

+ \\ii{t,x)+i$(t,x)\\ L p , (t,x) e [o,r] x R f/ . 

By a Gronwall-type estimate, which is possible because of A5 (see the proof of [5, 
Thm. 4.7(3)] for an elaboration of an argument of this type), we conclude 

sup v N (t,x)<C(T) sup \\Ii (t,x) + Ij (t,x)\\ LP . 

(t^)6[o,r]x/r (f,Jt)e[o,r]xA: 

where C(T) does not depend on K because of (11). For/^^x) we have for p > 1 

PfM \\w< 


+ 


< CI 1 + sup 

\ (f,x)e[o,r]x 


[ [ \G(t,x;s,y)-G N (t,x;s,y)](j(Y(s,y))B l (ds,dy) 

Jo Jr 1 ' 

f f \G(t,x;s,y)-G N (t,x;s,y)}(j{Y(s,y))M{ds,dy) 

Jo m d 

||F(f,x)||L/^ ^ J^\G{t,x;s,y)\\Bi\{ds,dy) 

+ ( f f |G(f,x;s,y)z| /, v(dj,dy,dz) N ) 

\Jo J(u N ) c J 


<C(T)r^(T,K) 


(26) 


uniformly in (f,x) £ [0. T] x K. In similar fashion one proves the estimate (26) for 
p g (0,1], perhaps with a different C(T). Next, when p > 1, (22) implies 
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KMII^= f [J, N G N (t,x;s,y)a(Y N (s,y))z(n-v)(d S ,dy,dz) 
Jo JR d J[-e N ,e N ] 


LP 


< 


I ^ ^A\Yy)o{Y N Y,y))zt {A ^ y} v{dsA yi dz) 

<40 ( ( [ [ , f N N |G(f,x;s,y)z| p 7r(d^)d( s ,y) N ) 

\ \J 0 J K rf J[-e N ,E N ] ) 

I / K J G(? ’ w)| Vm d 4)0j 


LP 


<C(T)(r» + r»). 


The case p £ (0,1] can be treated similarly, cf. the estimation of (t,x) above. 

It remains to prove that for each (t,x) £ [0. T] x K the convergence of Y N (t,x) 
to Y(t,x) is almost sure when r^(T,K), and are p*-summable. To this end, 
choose an arbitrary sequence (ajv)iVeN Q (0,1) converging to 0 such that 


Y An < 00 

N=l 


with An := 


(i#(T,K) + ri[ + r%y 


N 


Such a sequence always exists, see [13, Thm. 175.4], for example. So by (21) and 
Chebyshev’s inequality we derive 


P [|T(f,x) — Y N (t,x)\ > on] < 


\\Y(t : x)~Y N (t,x)\\ 

4 


p* 

LP 


< C(T)A N . 


Our assertion now follows from the Borel-Cantelli lemma. □ 


Example 2. The rates Y/ and from (17) only depend on the underlying Levy 
measure n. Let p,q £ (0,2] with q < p, and assume that /j_, jj \z\ q n[dz) < If 
A £ Yq, assume that q < 1. Then 


= 0([e N y p - q V p *^ , 





<0((e N y- q t {p>hAm ) . 


For instance, if e N = \/N k , then the sequence (r^Y* = &{N^ p ~ q ' 1 ) is summable 
for all k > (p — q)~ l . So in order to obtain a.s. convergence of Y N (t,x) —> Y(t ,x), 
a sufficient condition is to choose the truncation rates e N small enough. Similar 
conclusions are valid for the other two rates r^(T,K) and iY . 
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4 Truncation via Series Representations 

From the viewpoint of simulation, the truncation of the small jumps as presented in 
the previous section, has two main drawbacks: first, it may not be so easy to simulate 
the jumps of the truncated Levy measure, i.e. from the distribution 7r/7 t([— e,e] c ), 
for a large number of times; second, the jumps have to be simulated all over again 
when one goes from step N to step N + 1. These two problems can be overcome 
by using series representations for the Levy basis. The idea, going back to [16, 17] 
and already applied to the simulation of Levy processes [18], is to choose the jumps 
to be simulated in a random order. Instead of selecting the big jumps first and the 
smaller jumps later as in Sect. 3, we only choose the big jumps first more likely. 
The details are as follows: we fix a finite time horizon T £ R + and, recalling A7, a 
partition (Q')i'gN of R rf into pairwise disjoint compact sets such that U N = |J;Li Q‘■ 
We now assume that the jump measure p of A on the strip [0,7’] x R d x R can be 
represented in the form 


p(dt, dx. dz) = ^p,(dr,dx,dz), 

!= 1 

Hi(dt,dx,dz) = £ S (T i£ij 1 ( r iyi))(dt,dx,dz) a.s., (27) 

7=1 1 1 ' 1 

where H: (0,°°) x R —> R is a measurable function, satisfying H(-,v) = —//(•, — v) 
for all v £ R when A £ .9, and the random variables involved have the following 
properties for each i £ N: 

• (z l j: j £ N) and (^j: j £ N) are i.i.d. sequences with uniform distribution on 
[0,T] and Q', respectively. 

• (TV: j £ N) is a random walk whose increments are exponentially distributed 
with mean 1/77 

• (Vj: j £ N) is an i.i.d. sequence with distribution F on R, which we should be 

able to simulate from. We assume that F is symmetric when A £ . 

• The sequences z‘, c', f and V' are independent from each other. 

• (z , ,£\r , ,V l ) is independent from (z k ,% k ,r k ,V k : 

Because of (6), jX can always be written in the form (27) whenever the underlying 
stochastic basis is rich enough. We give three examples of such series representa¬ 
tions. 

Example 3. The proofs that the following choices are valid can be found in [18, 
Sect. 3], where also more examples are discussed. 

1. LePage’s method', we set F \= (5_i + 8\)/2 and H(r,±l) := ±p <_ (r,±l), 
where p^(r, ±1) = inf{x £ (0,°°): 7r(±[x,°o)) < r} for r £ (0,°°). 

2. Bondesson’s method', we assume that 7t(A) = f^°F(A/g(t))dt for A £ 3§(M. d ) 
with some non-increasing g: R + —» R + . Then we define H(r,v) := g(r)v. 

3. Thinning method: we choose F in such a way that Q is absolutely continuous 
with respect to F with density q, and define H(r,v) := V l{g( v )>r}- 
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Our approximation scheme is basically the same as in Sect. 3: we define G N by 
(16) and Y N by (19), with the difference that A N on [0,7’] x W ! is now defined as 

A A '(df,dr) := [ zjl N (dt,dx,dz) , 

J R 

/^(dLd^dz) := Y, E 5 (T':|'/r(r/v!))( dr,dr,ck ) ■ (28) 

i= 1 j:r[<N r r 1 


We can therefore rewrite Y N (t,x) for (t,x) E [0 ,T] x R^ as 


N 


Y N (t,x) = Y 0 (t,x) + '£ E G(t,x;T'^j)cr(F iV (Tj^j))77(77,y;)l {Tj<r} . (29) 
' 1 y: rf<N 


This yields the following simulation algorithm: 

Algorithm 2. Let ^ be a finite grid in [0, T] x R f/ and iVGN. 

1. For each 7 — 1,,7V set j := 1 and repeat the following: 

a. Draw £' from an exponential distribution with mean l/T. 

b. Define rj := +E i j (rf := 0). 

c. If rj > N, set Jj := j — 1 and leave the loop; otherwise set j := j + 1. 

2. For each i= 1 and j = 1..... J, simulate independently 

a. a pair (rj, ^j) with uniform distribution on [0,7’] x Q'\ 

b. a random variable Vj with distribution F\ 

c. the random variable 

3. Sort the sequence (rj: i = 1,... ,N,j = 1,... ,7,j in increasing order, yielding 

sequences (t,-, J], V,-: / ; ). Now define 


-l 


y^T,^) := T 0 (T^) + E T/, Zj)o(Y N {Tj,Zj))H(rj,Vj) . 

7=1 


4. For each (t,x) £ Sf simulate Yo(t,x) and define T(f ,x) by (29). 

We can now prove a convergence theorem for K v to T, similar to Thm. 1. Define 




(jT|ff(r,v)|'F(dv)dr) ^ , 

L X h ^ ,v ) iL{p>i ’ a ^} f ( dv ) dr 


(30) 


Theorem 2. Grant assumptions A1-A7 under which the SVE (1) has a unique solu¬ 
tion in /i[[ )c . Further suppose that the jump measure LI of A has a representation in 
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form of (27). Then Y N as defined in (29) belongs to B p oc for all N £N, and for all 
T £ R+ and compact sets K C R c/ there exists a constant C{T) £ R+ independent 
ofN and K such that 

sup \\Y{t : x)-Y N (t,x)\\Lr <C(T)(r J ?(T,K) + ri[ + r%) . (31) 

(f,jc)6[0,r]x/f 


IfL: =1 (r»(T,K) + r» + r$)P* <~, then we also have for all (t,x) £ [0, T] x K that 
Y N (t,x) —>• Y (; t,x ) a.s. as N —> °° 

Proof We start with some preliminaries. It follows from (27) and [18, Prop. 2.1] 
that on [0,r] x R f/ x R we have v = vo h~ l where v(dr,d.r,dr,dv) = df cbcdrF(dv) 
and h(t,x,r,v) = (t,x,H(r, v)). Therefore, conditions (6) and (9) imply that 

f f \H(r,v)\ p F(dv)dr= [ \z\ p 7l(dz) <°° , 

Jo Jr Jr 

and \H{r,v)\l {p>1A< £y } F{dv)dr= j^\z\t {p>XMy] n(dz) < °° . 

Consequently, rf and are well defined and converge to 0 when IV —>• Similarly, 
the compensator Vn of the measure p — p N is given by y.vfdf .tl.v.dz) = dtdxn A (dz), 
where n N = (Leb®F) oFJf 1 and FJ N (r,v) = H(r,v )1 (#,<*>) M- 

For the actual proof of Thm. 2 one can basically follow the proof of Thm. 1. 
Only the estimation of I^(t,x) as defined in (23) is different, which we shall carry 
out now. In the case of p > 1, we again use the Burkholder-Davis-Gundy inequality 
and obtain for (t.x) £ [0, F] xl rf 


\\I^{t,x)\\LP = [ j [ G N (t,x;s,y)o(Y N (s,y))z(p N -V N )(ds,dy,dz) 

JO JR rf JR 

+ [ [ [ G N {t,x;s,y)a(Y N {s,y))zt{ A ^ } V N (ds,dy,dz) 
JO JJ ]R 

<C(T)l ( ( [ [ [ \G(t,x-,s,y)H(r,v)\ p F(dv)drd{s,y) 

\ \J o Jm d Jn Jr 

+ [ / [ [ H (T v ) t {A^}F{dv)dr 

Jo JmJ Jn Jr ^ J 

<C(F)(/f + ^). 


LP 

l /P 


d (s,y) 


The case p £ (0,1] is treated analogously. One only needs to replace /Xy — Vy by tty 
and estimate via Jensen’s inequality. □ 

Example 4 (Continuation of Ex. 3). We calculate the rates r% and from (30) for 
the series representations given in Ex. 3. We assume that p.q £ (0,2] with q < p are 
chosen such that Jj_j jj |j| £/ 7r(dz) < °°, and q < 1 if A £ fig. For all three examples 
we use the fact that K = (Leb®F)o// _1 and that r > IV implies \H(r, v)| < |//(7V, v) | 
for all v £ R. 
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1. LePage’s method'. We have 

\H(r,l)\P+\H(r,-l)\P 


= r 

JN 


dr < 


d 




|z| P 7T(dz) 


< 


- [ 

2 J\H(l 


[H(l,-l),ff(l,l)] 


|z| 9 7r(dz)(|//(lV, —1)| V \H(N , l)|) p ~‘ ? , 


and therefore 


2. Bondesson’s method : Since II(r. v) = g(r)v and g is non-increasing, we obtain 


('f) P * = T [ \g(r)v\PF(dv)dr<(g(N)y^ f g*(r)dr [ |v| p F(dv) , 
jn Jr Jo Jr 

and consequently 


A =^(g(A0( p -«)/ p *) , r» = ff{g{N) X -n {p>lMy} ) . 

3. Thinning method'. Here we have 

fq{v)VN __ ,\ 1/p * / r . , n (q(v)-N)V0 ^ 1/p ‘ 


rd= 


< 




i Ip* 

i /p* 


r rq[v)\n\ 

/ / |v| p drE(dv) 

Jr Jn 

j R \z\ P H q { V )>N} n ^ 

^ j k|l{g(v)>iV} 1 {p>l,A^} ?I: (dz) ■ 


(f /’ y ’ " ?r(dv) 

VJr <?(v) 


In most situations, there exist (c^men Q with s N —» 0 as IV — > °° such that 

{tfO) > N} C [— e N ,e N ], In this case, one can apply the estimates in Ex. 2. 


5 Simulation Study 

In this section we visualize the sample path behaviour of the stochastic heat equation 
from Ex. 2 via a simulation study, using MATLAB programs from [4]. We take A 
to be a Levy basis without drift, whose Levy measure 7Z is that of a gamma process, 

i.e. 

n(dz) = yz _1 exp (-Az)l{ z>0 } dz 

with two parameters y,A>0. In the figures below their values are always y= 10 
and A = 0.1. Furthermore, we set Fq = 0 and <7=1. Especially the latter choice 
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N=50 


N=100 



Fig. 1 Successive approximations of Y as given in (32) via Bondesson’s method in dimension 1 
for ( t,x ) 6 [0,1] x [0,1] with N 6 {50,100,250,500} jumps in the region [0,1] x [—1,2] 


simplifies the subsequent discussion a lot, but none of the issues we address below 
relies on this assumption. Thus, the process we would like to simulate is 

Y (t,x) = f f g(t-s,x-y)A(ds,dy), (t,x) G R+ x R f ', (32) 

J0 J R rf 

with g being the heat kernel given in (14). In order to understand the path properties 
of Y, it is important to notice that g is smooth on the whole IB. | x W 1 except at the 
origin where it explodes. More precisely, for every t £ (0,°°) the function xi-tg(t,x) 
is the Gaussian density with mean 0 and variance 2 1, which is smooth and assumes 
its maximum at 0. Also, for every x ^ 0, the function t g(t,x ) is smooth (also at 
t = 0), with maximum at t = ||x|| 2 / (2d). However, if x = 0, then g(t, 0) = (4 nt)~ d ! 2 
has a singularity at 1 = 0. 

These analytical properties have direct consequences on the sample paths of Y. 
When A is of compound Poisson type, that is, has only finitely many atoms on 
compact sets, it can be readily seen from (32) that the evolution of Y after a jump J at 
(t, 4) follows the shape of the heat kernel until a next jump arrives. In particular, for 
x = %,Y(t,x) jumps to infinity at T, and decays in t like J(An(t — t )) - ' 7 / 2 afterwards. 
But for every x ^ £, the evolution 1 1 -> Y(t ,x) is smooth at t = T. In fact, it first starts 
to increase until t = t + ||x — E, || 2 / (2d) and then decays again. As a consequence, 
in space dimension 1, the space-time plot of Y shows a basically smoothly evolving 
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50 
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50 


0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 


Fig. 2 Several f-sections of the realization of Y shown in Fig. 1 with N = 500 


path, only interrupted with slim poles at the jump locations of A ; see the case N = 50 
in Fig. 1. However, when A has infinite activity, that is, has infinitely many jumps 
on any non-empty open set, then it is known from [16, Thm. 4] that on any such set 
Y is unbounded, at least with positive probability. Therefore, the space-time plots 
of the approximations of Y with finitely many jumps must be treated with caution: 
in the limiting situation, no smooth area exists any more, but there will be a dense 
subset of singularities on the plane, which is in line with Fig. 1. 

Another interesting observation, however, is the following: if we consider a 
countable number of x- or f-sections of Y (forx £ M. d , the x-section of Y is given by 
the function f 1 —s- Y (f ,x); for f £ K + , the f-section of Y is the function n-tf (f ,x)), 
then it is shown in [19, Sect. 2] that these are continuous with probability one. In¬ 
tuitively, this is possible because a.s. the sections never hit a jump (although they 
are arbitrarily close). For instance. Figs. 2 and 3 show f-sections of a realization of 
(32) in one, respectively two space dimensions. So as long as we only take count¬ 
ably many “measurements”, we do not observe the space-time singularities of Y 
but only its relatively regular sections. In theory, this also includes the x-sections of 
the process Y. But if we plot them for one space dimension as in Fig. 4, one would 
conjecture from the simulation that they exhibit jumps in time. However, this is not 
true: the jump-like appearance of the x-sections are due to the fact that g(-,x) resem¬ 
bles a discontinuous function at f = 0 for small x. Of course, it follows right from 
the definition (14) that all x-sections of g are smooth everywhere. 
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t=0.28 t=0.29 1=0.30 



Fig. 3 Several /-sections in the region [— 1 , l] 2 of a realization of Y in dimension 2 by Bondesson’s 
method with N = 500 jumps within [0,1] x [—2,2] 2 




Fig. 4 The x-section of the realization of Y as in Fig. 1 with N = 500 at x = 0.6 and the heat kernel 
g(-,x) at x = 0.002 
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6 Conclusion and Outlook 

In Sects. 3 and 4 we have presented two simulation algorithms for the SVE (1): 
Algorithms 1 and 2. In Thms. 1 and 2 we have determined the rate of convergence of 
the approximations Y N to Y in the IP- sense. If these rates are small enough, we have 
also proved a.s. convergence. Although the theoretical analysis of both schemes lead 
to quite similar results regarding their convergence behaviour, there are important 
differences which will decide on whether the one or the other method is preferable in 
concrete situations. For the first method of truncating the small jumps to work, one 
must be able to efficiently simulate from the truncated Levy measure n/n(\— £,e] c ) 
for small e. For the second method, which relies on series representations, the main 
challenge is to choose H and F in a way such that H is explicitly known and F can 
be easily simulated from. For instance, if one uses LePage’s method (see Ex. 3), then 
F = (S | + 5i)/2 is easily simulated, but for H, which is given by the generalized 
inverse tails of the underlying Levy measure, maybe no tractable expression exists. 

Finally, let us comment on further generalizations of the our results. Throughout 
this paper, we have assumed that the driving noise A is a homogeneous Levy ba¬ 
sis, i.e. satisfies (5). In fact, we have introduced this condition only for the sake of 
simplicity: with a straightforward adjustment, all results obtained in this paper also 
hold for time- and space-varying (but deterministic) characteristics. Another issue 
is the finite time perspective which we have taken up for our analysis. An interest¬ 
ing question would be under which conditions (1) has a stationary solution, and in 
this case, whether one can simulate from it. Sufficient conditions for the existence 
and uniqueness of stationary solutions to (1) are determined in [5, Thm. 4.8]. Under 
these conditions, the methods used to derive Thms. 1 and 2 can indeed be extended 
to the case of infinite time horizon. We leave the details to the reader at this point. 

At last, also the hypothesis that A is of pure-jump type can be weakened. If A 
has an additional drift (including the case where A has locally infinite variation and 
is not symmetric) but still no Gaussian part, the approximations Y N in (19) or (29) 
will contain a further term that is a Volterra integral with respect to the Lebesgue 
measure. So each time in between two simulated jumps, a deterministic Volterra 
equation has to be solved numerically, which boils down to a deterministic PDE 
in the case where G comes from an SPDE. For this subject, there exists a huge 
literature, which is, of course, also very different to the stochastic case as consid¬ 
ered above. If A also contains a Gaussian part, then one has to apply techniques 
from the papers cited in Sect. 1 and ours simultaneously. We content ourselves with 
referring to [22], who numerically analyzes a Volterra equation driven by a drift 
plus a Brownian motion. Finally, let us remark that if p = 2 (in particular, G must 
be square-integrable), it is possible for some Levy bases to improve the results of 
Sect. 3 if we do not neglect the small jumps completely but approximate them via a 
Gaussian noise with the same variance, cf. [1] in the case of Levy processes. 
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